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Abstract 

A parallel dispersive finite-difference time-domain (FDTD) method for the 
modeling of three-dimensional (3-D) electromagnetic cloaking structures is 
presented in this paper. The permittivity and permeability of the cloak are 
mapped to the Drude dispersion model and taken into account in FDTD sim- 
ulations using an auxiliary differential equation (ADE) method. It is shown 
that the correction of numerical material parameters and the slow switching- 
on of source are necessary to ensure stable and convergent single-frequency 
simulations. Numerical results from wideband simulations demonstrate that 
waves passing through a three-dimensional cloak experience considerable de- 
lay comparing with the free space propagations, as well as pulse broadening 
and blue-shift effects. 
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1. Introduction 

Recently, a great deal of attention has been paid to the analysis and de- 
sign of electromagnetic cloaking structures, since first proposed by Pendry 
et al. [H. The specially designed cloak is able to guide waves to propa- 
gate around its central region, rendering the objects placed inside invisible 
to external electromagnetic radiations. The ideal cloak in pLj] requires inho- 
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mogeneous and anisotropic media, with both permittivity and permeabihty 
independently controlled and radially dependent, making its practical real- 
ization very difficult. Therefore it has been proposed to use simplified mate- 
rial parameters for both transverse electric (TE) ^ and transverse magnetic 
(TM) g cases. In order to reduce the scattering due to the impedance 
mismatch introduced by the simplified cloaks, an improved linear cloak a 
high-order transformation based cloak j^, and a 'square root' transformation 
based cloak have also been proposed. 

The coordinate transformation technique used in [l|, 0] has also been 
applied to the design of magnifying perfect and super lenses jsf, electro- 
magnetic field rotators HL the reflectionless complex media for shifting and 
splitting optical beams [lO|, and conformal antennas ll| etc. The experi- 



mental demonstration of a simplified cloak consistin g of split-ring resonators 



(SRRs) has been reported at microwave frequencies [12|. For the optical fre- 
quency range, the cloak can be constructed by embedding silver wires in a 
dielectric medium jsf, or using a gold-dielectric concentric layered structure 
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The modeling of Pendry's invisible cloak has been performed by using 
both analytical and numerical methods. Besides the widely used coordi- 
nate transformation technique [H, Q, H, [o], 0, 11, 15 1, a cylindrical wave 
expansion technique 
tering model 
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and a method based on the full-wave Mie scat- 
171 . |18| have also been applied. In addition, the full-wave 
finite element method (FEM) based commercial simulation software COM- 
SOL Multiphysics"""^ has been extensivelyused to model different cloaks and 
validate theoretical predictions B a a 

y, ll9| , due to its ability of dealing 
with anisotropic and radial dependent material parameters. However similar 
to other frequency domain techniques, the FEM may become inefficient when 
wideband solutions are needed. So far, the time domain techniques that have 
been developed to model the cloaking structures include the time-dependent 
scattering theory 20], the transmission line method (TLM) j2l| and the 
finite-difference time domain (FDTD) method 22|. Due to its simplicity in 



implementation and the ability of treating anisotropic, inhomogeneous and 
nonlinear materials, the FDTD method has been extremely popular for the 
analysis of electromagnetic structures. However due to the computational 
complexity, so far the FDTD modeling of three-dimensional (3-D) cloaking 
structures has not been attempted. In this paper, we extend our previously 
proposed 2-D FDTD method [22*] to the 3-D case and develop a parallel 
dispersive FDTD method to model 3-D cloaking structures and reveal the 
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extraordinary behavior that is different from its counterparts under 2-D as- 
sumptions. 



2. Parallel dispersive FDTD modeling of 3-D cloaks 

A complete set of material parameters of the ideal cloak in spherical 
coordinate is given by [l| 



R2 (r-Ri^^ 



R2 ~ Ri 

R2 
R2 ~ Ri 



_ _ R2 

K2 — rCi 

where Ri and R2 are the inner and outer radii of the cloak, respectively, 
and r is the distance from a spatial point within the cloak to the center of 
the cloak. Since the values of Er and fir are less than one (between and 
{R2 — Ri) /R2), same as the case for left-handed materials (LHMs), the cloak 
cannot be directly modeled using the conventional FDTD method. However, 
one can map the material parameters using dispersive material models, for 
example, the Drude model 

^2 

Briu) = 1 - , r , (2) 

where Up and 7 are the plasma and collision frequencies of the material, 
respectively. By varying the plasma frequency, the radial dependent mate- 
rial parameters in ([T]) can be achieved. For example, for the ideal lossless 
case considered in this paper, i.e. the collision frequency in ([2]) is equal to 
zero (7 = 0), the radial dependent plasma frequency can be calculated using 
Up = — Er with a given value of Er calculated from ([1]). Note that in 
practice, the plasma frequency of the material depends on the periodicity 
of the split ring resonators (SRRs) or wires [3|, which varies along the 
radial direction. It should be noted that other dispersive material models 
(e.g. Debye, Lorentz etc.) can be also considered for the modeling of elec- 
tromagnetic cloaks. However, the Drude model has the simplest form when 
implemented using the dispersive FDTD method and has been widely used 
in the modeling of metamaterials. The Lorentz model can be also used in 
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FDTD simulations with some additional modifications to the iterative equa- 
tions. The Debye model may be less accurate to characterize the dispersion 
behavior of metamaterials and rarely used in the community. 

Since the conventional FDTD method 23|, |2J] deals with frequency- 



independent materials, the frequency-dependent FDTD method is hence re- 
ferred as the dispersive FDTD method 25|, l26|, l27|]. There are also differ- 



ent dispersive FDTD methods using different approaches to deal with the 
freq uency-dependent material parameters: the recursive convolution (RC) 
method j25|, the auxiliary differential equation (ADE) method 26[] and the 



Z-transform method [27j. Due to its simplicity, we have chosen the ADE 
method to model 3-D cloaks in this paper. 

The ADE dispersive FDTD method is based on the Faraday's and Am- 
pere's Laws 

(3) 



V X E 



V X H 



dt 



(4) 



eE and B = /iH where e and jj, 
(jlj) can be discretized following a 
standard procedure [23|, l2j] , which leads to the conventional FDTD updating 
equations 



as well as the constitutive relations D = 
are expressed by ([1]). Equations ([3]) and 



At • V X E"+5, 
At- V X H"+i 



(5) 
(6) 



where V is the discrete curl operator. At is the discrete FDTD time step and 
n is the number of the time steps. 

In addition, auxiliary differential equations need to be taken into account 
and they can be discretized through the following steps. For the conventional 
Cartesian FDTD mesh, since the material parameters given in are in 



spherical coordinates, the following coordinate transformation is used [28 
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The tensor form of the constitutive relation is given by 
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Note that the inverse of the permittivity tensor matrix ([9]) exists only when 
7^ 0, £0 7^ and 7^ 0. However the inner boundary of the cloak does not 
satisfy the condition of 7^ 0. Therefore in our FDTD simulations, we place 
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a perfect electric conductor (PEC) sphere with radius equal to Ri inside the 
cloak to guarantee the validity of (jH]). 
Substituting Qj into ([8]) gives 
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Since the above equations have a similar form, in the following, the derivation 
of the updating equation is only given for the component. The updating 
equations for the By and components can be derived following the same 
procedure. 

Express Er in the Drude form of ([T]), Eq. flTT]) can be written as 
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Note that ee and e,^ remain in (HM because their values are always greater 
than one and can be directly used in the conventional FDTD updating equa- 
tions 23, 2^. Applying the inverse Fourier transform and the following rules: 
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Eq. (1141) can be rewritten in the time domain as 
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The FDTD simulation domain is represented by an equally spaced 3-D 
grid with the periods Ax, Ay and Az along the x-, y- and z-directions, 
respectively. For the discretization of Eq. (fT6|) . we use the central finite 
difference operators in time {6t and 6^) and the central average operator 
with respect to time {fit and /x^): 



d 6t 
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where the operators 6t, S'f, Ht and /i^ are defined as in [29 
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where F represents field components and m^jmy^rriz are the indices corre- 
sponding to a certain discretization point in the FDTD domain. The dis- 
cretized Eq. f|T6|) reads 
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Note that in (fTSll . the discretization of the term Up of ( fT6l) is performed using 
the central average operator fif in order to guarantee the improved stability; 
the central average operator fit is used for the term containing 7 to preserve 
the second-order feature of the equation. Equation f|T8|) can be written as 
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After simple manipulations, the updating equation for can be obtained 
as 
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+b^xzDT^' + + - (ai,£;: + a2.£;r') l/«o., (20) 



where the coefficients are given by 
1 



o-Ox — £0 



<^2x = £0 



'JOxy 



1 



7 




2At 


4 






7 





2 , 
(At)2 ^ 2 



^OxK = sin^ cos^ ( 



sin^ 9 cos^ 



(At)2 ^ 2At 

, 2 



+ 



cos^ 6* cos^ (j) ^ siiP I 



b2xx = sin^ cos^ ( 



(At)2 

1 



+ 



cos^ cos^ 6 sin^ ( 



+ 



S4, 



2 , 
(Ai)2 ^ 2 



1 7 
(At)2 2At T 

2 



sin^ 9 sin cos 



2_ 

(At)2 2At 
1 



cos2 9 cos2 ^ sin^ ( 



£4, 



_J: 2L + i^ 

(Ai)2 2At 4 



(At)2 2At 



cos ^ sin 6 cos c!> sin c!) cos < 



S0 



£4, 



7 , 



(At)2 ' 2At 4 



bixy = — Sin" y sin cj) cos < 



^2x3/ = sin 0sin^cos( 



(At)2 

1 



+ 



cos2 sin (p cos </> sin cos (p 



£0 



(Ai)^ 



(At)2 2Ai 



+ 



cos 9 sin 6 cos (A sin S cos < 



£e 



1 1 . 



(Ai)2 2 At 4 
1 



boxz = sin ^ cos ^ cos 



(At)2 2At 



sin 9 cos ^ cos ( 

£9 



£<^ 



i 7 
(At)2 2At T 



9 



bixz = — sin ^ cos 6 cos < 



sin y cos y COS ( 



(Ai)2 

1 



Sin y COS tf cos t 



2At 



sin cos cos < 



^ + ^ 

(At)2 2At 4 



Following the same procedure, the updating equation for Ey is 
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Note that the field quantities Dr^, Dy and in are locally aver- 

aged values of D^^, Dy and D^, respectively since the x-, y- and 2;-components 
of the electric fields are in different locations in the FDTD domain 30|. How- 
ever, the averaging needs to be applied along different directions depending 
on the updating equations. Specifically in (l20l) . the averaged Dy and Dz can 
be calculated using 

— Dy{l,J,k)+Dy{t+l,J,k)+ Dyit,J -l,k)+Dyit + l,J -l,k) 

Uy[l,J,K) — 



Dzii,j,k) 



4 

Dz{i, i, k)+Dzit + 1, J, k)+Dzii, j,k-l)+Dzit + l,j,k~l] 



where k) is the coordinate of the field component. In fl2Tl) . the averaged 
D^ and D^ can be calculated using 

D,{i,j,k) + D,{i,j + l,k) + D,{i-l,j,k) + D,{i-l,j + l,k) 
D^{i,],k) = , 

_ . _ Dz{i,j,k) + Dz{i,j + l,k) + Dz{i,j,k-l)+Dzit,j + l,k-l) 

And in fl22|) . the averaged D^ and Dy can be calculated using 

D,(z,j, + fc + l) + D,(z-l,j,A;) + Z},(z-l,j, A; + l) 
Dx{t,j,k) = , 

— . _ Dy{Z,J, k)+Dy{t,j,k + l)+Dy{l,J -l,k)+Dy{l,J -l,k+l) 
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The updating equations for the magnetic fields if^., Hy and are in 
the same form as fl2Ul) - fl2^ with the same coefficients, and can be obtained 
by replacing E with H and D with B. The averaged field components can 
be calculated in a similar manner. Equations ([5]), ([6]), (!20l) - (!22l) and the 
updating equations for H from B (not given) form the updating equation set 
for the modeling of 3-D cloaks using the well-known leap-frog scheme 23|]. If 
the plasma frequency in ([21) is equal to zero i.e. ojp = 0, and se = jJ^e = £<t> = 

= 1, the above updating equation set reduces to the updating equation 
set for the free space. 

Since the FDTD method is inherently a numerical technique, the spatial 
as well as time discretization has important effects on the accuracy of simu- 
lation results. Also because the permittivity and permeability are frequency 
dependent, one can expect a slight difference between the analytical and 
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numerical material parameters due to the discrete time step in the FDTD 



method. From our previous analysis jSli that for the modeling of metamate 



rials especially the LHMs, the numerical errors due to the time discretization 
will cause spurious resonances, hence a req uirement of Ax < A/80 is neces- 



sary. Following the same approach as in [2^ , one can find that the numerical 



permittivity for the ideal 3-D cloak takes the following form: 

(23) 
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Note that Eq. fl23|) simplifies to the analytical Drude dispersion model ([2]) 
when At 0. With the expression of the numerical permittivity fl2^ avail- 
able, one can correct the errors introduced by the discrepancy between the 
numerical and analytical material parameters. For example, if the required 
permittivity is = e[. + je", after simple derivations, the corrected plasma 
and collision frequencies can be calculated as 
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The above averaging of field components and the correction of numerical 
material parameters ensure stable and accurate FDTD simulations of the 3-D 
cloak. However if the averaging is not applied, the field distribution becomes 
unsymmetrical and hence incorrect; if the numerical material parameters 
are not corrected, the FDTD simulations become unstable before reaching 
the steady-state. Therefore in our simulations of the 3-D cloaks, the field 
averaging and corrected material parameters are always used. 

The FDTD method is a versatile numerical technique. However, similar 
to other numerical methods, it is computationally intensive. For large elec- 
tromagnetic problems such as the modeling of 3-D cloaks, the requirement 
for system resources is beyond the capability of a single personal computer 
(PC). One way to resolve this problem is to divide the whole computational 
domain into many smaller sub-domains, and each sub-domain can be handled 
by a PC. By linking the PCs together with an appropriate synchronization 
procedure, the original large problem can be decomposed and solved effi- 
ciently. 
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One of the most attractive features of the FDTD method is that it can 
be easily parallehzed with very httle modifications to the algorithm. Since 
it solves Maxwell equations in the time-space domain, the parallel FDTD 

The data 
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algorithm is based on the space decomposition technique [33|, 
transfer functionality between processors (PCs) is provided by the message 
passing interface (MPI) library. Data exchange is required only for the ad- 
jacent cells at the interface between different sub-domains and is performed 
at each time step, hence the parallel FDTD algorithm is a self-synchronized 
process. Figure [1] shows the arrangement of the field components in different 
sub-domains in parallel FDTD simulations. The red arrows are the trans- 



-> 



XEx ^^^^ 
Sub-domain I 

























Ey 









Sub-domain II 



Figure 1 : The arrangement of field components in different sub-domains in parallel FDTD 
simulations. The red arrows indicate the field components which are transferred from the 
neighboring domain during the data communication process and used to update the field 
components on the boundary of the current sub-domain. 



ferred field components from the neighboring sub-domain during the data 
communication process. At the end of parallel FDTD simulations, the re- 
sults calculated at each processor need to be combined to obtain the final 
simulation result. In comparison to conventional parallel FDTD method, 
the parallelization of the dispersive FDTD method introduced in this pa- 
per requires additional field components to be transferred between adjacent 
sub-domains during the synchronization process, because of the applied field 
averaging scheme. The complexity of algorithm further increases if the whole 
computational domain is divided along more than one direction, although the 
data communication load and the overall simulation time may be reduced. 

The PC cluster used to simulate 3-D cloaks in Queen Mary, University of 
London consists of one head node for monitoring purposes and fifteen com- 
pute nodes for performing calculation tasks. Each node has Dual Intel Xeon 
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E5405 (Quad Core 2.0 GHz) central processing units (CPUs) and there are 
128 cores and 512 GB memory in total. The nodes are connected by a 24- 
port gigabit switch. The GNU C compiler (GCC) and a free version of MPI, 



MPI Chameleon (MPICH), developed by Argonne National Laboratory 32 



are used to compile the developed parallel dispersive FDTD code and han- 
dle the inter-core data communications, respectively. The above developed 
parallel dispersive FDTD method has been implemented to model the ideal 
3-D cloaks, and the simulation results and discussions are presented in the 
following section. 



3. Numerical results and discussions 




Figure 2: The 3-D parallel dispersive FDTD simulation domain for the case of plane- wave 
incidence on the cloak. The red rectangle indicates the location of the source plane. 



size in all simulations is Ax = Ay = A/ 150 where A is the wavelength at 
the operating frequency / = 2.0 GHz. The time step is chosen according 
to the Courant stability criterion [[2^] i.e. At = Ax/y^c The radii of the 
cloak are i?i = 0.1 m and i?2 = 0.2 m. In the present paper, only the 
ideal case (lossless) is considered i.e. the colhsion frequency in ((21) is equal 
to zero (7 = 0). The radial dependent plasma frequency can be calculated 
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using with a given value of Sj. calculated from (JT]). The computational 
domain is truncated using Berenger's perfectly matched layer (PML) [s^ 
in ^/-direction to absorb waves leaving the computational domain without 
introducing reflections, and terminated with periodic boundary conditions 
(PBCs) in X- and z-directions for the modeling of a plane-wave source. The 
plane-wave source is implemented by specifying a complete plane of FDTD 
cells using a certain wave function. The electric and magnetic fields of the 
plane wave are along the z- and respectively and the propagation 

direction is along the y-axis, as indicated in Fig. [2l For simplicity, the whole 
simulation domain is only divided along ?/-direction into 100 sub-domains 
and in total 100 processors and 220 gigabyte (GB) memory were used to run 
the parallel dispersive FDTD simulations. Each simulation lasts around 45 
hours (13000 time steps) before reaching the steady-state. 

Figures [3] and H] show the normalized steady-state field distributions for 
the Ez and components in y-z and x-y planes, respectively. It can be seen 
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Figure 3: Normalized field distributions for the Ez component in (a)-(c) y-z plane and 
(d)-(f) x-y plane in the steady-state of the parallel dispersive FDTD simulations. The 
cutting planes are (see Fig. [2]): (a) x = 2A, (b) x = 4A/3, (c) a; = A, (d) z = 2A, (e) 
z = 4A/3, (f) z = X. The wave propagation direction is from left to right. 
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Figure 4: Normalized field distributions for the component in (a)-(c) y-z plane and 
(d)-(f) x-y plane in the steady-state of the parallel dispersive FDTD simulations. The 
cutting planes are (see Fig. [2]): (a) x — 2A, (b) x — 4A/3, (c) x = \, {d) z — 2A, (e) 
z — 4A/3, (f) z = \. The wave propagation direction is from left to right. 

that the plane wave is guided by the cloak to propagate around its central 
region, and recomposed back after leaving the cloak. There is nearly no 
reflection (except those tiny numerical ones due to the finite spatial resolution 
in FDTD simulations), since the material parameters ([1]) vary continuously 
in space while keeping the impedance the same as the free space one. It 
is also interesting to notice that the Ez component in y-z and x-y planes 
in Fig. [3] and the Hj. component in x-y and y-z planes in Fig. H] have the 
same distributions (with different amplitude), which is due to the fact that 
the ideal 3-D cloak is a rotationally symmetric structure with respect to the 
electric and magnetic fields. The wave behavior near the 3-D cloak can be 
better illustrated using the power flow diagram, as plotted in Fig. O It is 
shown that the Poynting vectors are diverted around the central area enclosed 
by the cloak. Therefore objects placed inside the cloak do not introduce any 
scattering to external radiations and hence become 'invisible'. 

The above presented results validate the developed parallel dispersive 
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Figure 5: Power flow diagram of a plane wave incidence on the ideal 3-D cloak calculated 
from parallel dispersive FDTD simulations. 

FDTD method and demonstrate the cloaking property of the structure. How- 
ever, there are some numerical issues that need to be addressed in FDTD sim- 
ulations. Besides the correction of numerical material parameters introduced 
earlier, since the cloak is a sensitive structure, for single-frequency simula- 
tions, the switching time of the sinusoidal source also has significant impact 
on the convergence time. Normalized field distributions from the simulations 
using different switching time are plotted at the time step t = 1320At and 
shown in Fig. [6l It can be seen that if the source is switched to its maximum 
amplitude within a short period of time, because of the multiple frequency 
components excited, and the cloak is essentially a narrowband structure due 
to its dispersive nature, the scattering from the cloak may occur, as shown in 
Fig. ^a). The scattered waves oscillate within the lossless cloak and hence 
it requires a very long time for the simulations to reach the steady-state. It 
is also demonstrated that if the switching time is greater than IOTq where Tq 
is the period of the sinusoidal wave, the scattered waves can be significantly 
reduced and a much shorter convergence time in simulations can be achieved. 
Therefore in the previous simulations, the switching time of SOTq was used. 

Since the FDTD method is a time domain technique, it is convenient to 
study the transient response of the 3-D cloak. The snapshots of the field 
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Figure 6: Comparison of the influence of different switching time (ST) of the sinusoidal 
source on the simulation results: (a) ST = Tq, (b) ST = IOTq, (c) ST = SOTq, where Tq 
is the period of the sinusoidal wave. The wave propagation direction is from left to right 
and the normalized field distributions are plotted in the x-y plane (z — 2A, see Fig. [5]) and 
at the time step t = 1320Ai. 

distributions for the component at different time steps t = 3000 At (5.77 
ns), t = 5000 At (9.62 ns), t = 8000 At (15.40 ns) are taken and plotted in 
Fig. [71 It is shown in Fig. [7](a) that outside the shadow region behind the 
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Figure 7: Snapshots of the field distributions for the Ez component at different time steps 
in the parallel dispersive FDTD simulations: (a) t = SOOOAt (5.77 ns), (b) t = 5000At 
(9.62 ns), (c) t = SOOOAi (15.40 ns), plotted in the x-y plane {z = 2A, see Fig. [2]). The 
wave propagation direction is from left to right. 

cloak {y ~ 3.5A, x < 0.5A and x > 3.5A), waves propagate at the speed 
of light and the wave front remains the same as the one before reaching 
the cloak. However due to the fact that the waves that travel through the 
cloak undergo a longer path compared to the free space one, and since the 
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group velocity cannot exceed the speed of light, the wave front experiences a 
considerable delay in forming back to the free space one in the shadow region 
behind the cloak, as it is illustrated by the field distributions at different time 
steps in FDTD simulations in Fig. [71 In fact, the convergence of simulations 
is quite slow and the steady-state is reached in simulations at around 13000 
time steps (25.02 ns). 

Another advantage of the FDTD method is that a wideband frequency 
response can be obtained with a single run of simulations. In comparison to 
the previous single-frequency case, a wideband Gaussian pulse with the cen- 
tral frequency of 2.0 GHz and covering the frequency range of 1.5 ~ 2.5 GHz 
is used instead. At 5 FDTD cells away at both the front and back of the cloak 
along its central axis {x = z = 2A, see Fig. [2]), the amplitude of is recorded 
during the simulation and then transformed to the frequency domain. The 
recorded time domain signals and their spectra are plotted in Fig. [HI It is 




2 4 6 8 10 12 1 1.5 2 2.5 3 

Time (ns) Frequency (GHz) 

Figure 8: (a) The recorded time domain signals at 5 FDTD cells away at both the front 
and back of the cloak along its central axis {x — z — 2A, see Fig. [2|). (b) The spectra of 
the recorded time domain signals. 

found from Fig. [8](a) that the time delay between the directly received time 
domain signal in front of the cloak (red solid line) and the received time 
domain signal through the cloak (blue dashed line) is approximately 3.4 ns. 
However the distance between the two recording points is 0.402 m and the 
expected time delay for a plane wave propagating in the free space is 1.34 
ns. This clearly demonstrates that the 3-D cloak introduces a time delay to 
the waves propagating through it. It is also found from Fig. [Ht^a) that the 
width of the pulse passing through the cloak has been broadened, which is 
due to the dispersive nature of the cloaking material. It is interesting to note 
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that in Fig. MJo), the spectrum of the directly received time domain signal is 
centered at 2.0 GHz, however the central frequency of the signal that passes 
through the cloak is considerably higher (~ 2.1 GHz). This frequency shift 



has been demonstrated theoretically in [36|, which is explained as that the 
frequency components higher than the working frequency of the cloak are 
enhanced and the frequency components lower than the working frequency 
of the cloak are weakened. The shift is found to be much more pronounced 
from our analysis since the observation point is taken near the back surface 
of the cloak in our simulations, and it was 30A away from the cloak consid- 



ered in [36| . The difference may be also due to the Lorentz dispersion model 
considered for the permeability of the cloak in [36] , while in our simulations, 
both the permittivity and permeability are assumed to follow the same Drude 
dispersion model. 



4. Conclusion 

In conclusion, a parallel dispersive FDTD method has been developed to 
model the ideal 3-D cloak. The radial dependent permittivity and perme- 
ability of the cloak are mapped to the Drude dispersion model and taken 
into account in FDTD simulations using an ADE based method. Due to the 
memory restraint of a single PC, a parallel FDTD method is developed to 
handle the large amount of memory and simulation time required to model 
the 3-D cloak. FDTD simulation results are validated by those obtained 
using analytical methods. It is demonstrated that for single-frequency simu- 
lations, the source excitation needs to be switched on slowly enough to avoid 
the wave scattering from the cloak, due to the sensitivity of the cloaking ma- 
terial. It is also shown from the transient FDTD analysis that waves passing 
through the cloak experience a considerable time delay comparing with the 
free space propagations, as well as a pulse broadening effect and a blue-shift 
of the pulse's central frequency. 

The ideal 3-D spherical cloak is considered in this paper. The method 
developed here can be also used to model cylindrical cloaks and compare 
their properties with the spherical one, such as the direction dependency 
issue. The ideal cloaks are lossless and the cloaking material properties vary 
in space continuously, which make their practical realization very difficult. 
The developed FDTD method can be also applied to study the effect of losses 
in cloaks, evaluate the performance of simplified cloaks as well as assist the 
design and realization of practical cloaks. 
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